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An asymptotic method for finding instabilities of arbitrary d-dimensional large-amplitude pat- 
terns in a wide class of reaction-diffusion systems is presented. The complete stability analysis of 
\£) , 2- and 3-dimensional localized patterns is carried out. It is shown that in the considered class of 

■ systems the criteria for different types of instabilities are universal. The specific nonlinearities enter 

the criteria only via three numerical constants of order one. The performed analysis explains the 
self-organization scenarios observed in the recent experiments and numerical simulations of some 
' concrete reaction-diffusion systems. 
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i—l ■ I. INTRODUCTION 

> ■ 

In the last two decades the problem of pattern formation and self-organization has become a paradigm of modern 
science p]-p|]. Patterns are observed in a wide variety of physical systems, such as gas and electron-hole plasma; 
various semiconductor, superconductor and gas-discharge structures; some ferroelectric, magnetic, and optical media; 
systems with uniformly generated combustion material (see [^|-0J^] and references therein). Pattern formation and 
self-organization are most conspicuous in chemical and biological systems (see and references therein). As a rule, 
all these systems are extremely complicated. In order to describe pattern formation phenomena in them a number of 
simplifications is made. The majority of the simplified models reduce to a pair of reaction-diffusion equations of the 
c/3 ' activator-inhibitor type 

(— ' ■ 
(— > , 
Ctf . 

Or 



On 



X 



T e —=l 2 A9-q(9,r ] ,A), (1) 
at 



r ri ^=L 2 A v -Q(9, v ,A), (2) 

where 9 and rj are the distributions of the activator and the inhibitor, respectively; q(9, r\, A) and Q{9, rj, A) are certain 
nonlinear functions; I and L are the characteristic length scales, and Tg and are the characteristic time scales of 9 
and 77, respectively; A is the control parameter. The well-known models of certain autocatalytic reactions, such as the 
Brusselator Q| , the two-component version of the Oregonator and the Gray-Scott jl(J] models, the classical model 
of morphogenesis proposed by Gierer and Meinhardt pT] ; the piecewise-linear Jl^,[l3| and FitzHugh-Nagumo JTi 15 



models describing the propagation of impulses in the nerve tissue are the special cases of Eqs. (Jl|) and (||). These 
models are most widely used in the analytical investigations of different types of patterns |0-Eq| . 

The main self-organization phenomenon in the considered systems is spontaneous transformation of one type of 
pattern to the other as certain parameters of the system are varied. Self-organization scenarios become extremely 
diverse in 2- and 3-dimensional systems. In this situation the most interesting are the spontaneous transformations 
of simple patterns, especially localized steady patterns — autosolitons (AS), into much more complicated ones. A 
lot of different types of these transformations were recently observed in the experiments with semiconductor and gas 



discharge structures |27-29|, chemical reaction-diffusion systems ]30|-ffi|, and in the numerical simulations |33] p5[ . 
At the same time, the general theory of such transformations is absent. 

The general asymptotic method of constructing the solutions in the form of large-amplitude localized, periodic, and 
more complex one-dimensional patterns, and the qualitative analysis of their stability in the systems described by 
Eqs. (0) and @ were developed by Kerner and Osipov More formal and mathematically rigorous analysis 
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of one-dimensional patterns was carried out by Nishiura, Mimura, and the co-authors [|43|-|45| . Static, pulsating, and 
traveling large-amplitude patterns in simple reaction-diffusion systems were studied in detail E , E^ - glL^5 12 



In two and three dimensions Kerner and Osipov have constructed the asymptotic solutions for radially-symmetric 
patterns, and also analyzed the stability of one-dimensional patterns in higher dimensions p,p|,p9|-p2l ■ Ohta, Mimura, 
and Kobayashi have developed an approach which allowed them to study the stability of one-dimensional and radially- 
symmetric patterns in 2- and 3-dimensional versions of simple piecewise- linear model of reaction-diffusion system [[47| . 
This approach was further developed by Petrich and Goldstein, who applied it to a version of the FitzHugh-Nagumo 
model |23|| , 

In this paper we develop a systematic procedure of finding the bifurcation points of an arbitrary d-dimensional 
pattern. Using this procedure, we analyze the stability of the major types of patterns in arbitrary dimensions. On 
the basis of this analysis, we draw conclusions about possible scenarios of the transformations of patterns. 

Our paper is organized as follows. In Sec. II we generalize the method of constructing the asymptotic solution for 
one-dimensional and radially-symmetric patterns developed in Refs. [pPj3r^-^2]| to the case of arbitrary d-dimensional 
patterns. In Sec. Ill we present the derivation of the general dispersion relation governing the stability of an arbitrary 
pattern. In Sec. IV we apply the obtained results to the one-dimensional AS in higher dimensions. In Sec. V we 
analyze the instabilities of the spherically-symmetric AS in three dimensions, and cylindrically-symmetric AS in two 
and three dimensions. In Sec. VI we summarize the obtained results and discuss their implications on the evolution 
of patterns and also give comparisons with the experimental and numerical data. 

II. ASYMPTOTIC SOLUTIONS FOR ARBITRARY D-DIMENSIONAL PATTERNS 

If we choose L and t v as the units of length and time, we can write Eqs. ([[]) and (Q) as 

a^ = e 2 A9-q(6,r l> A), (3) 



^=Ar,-Q(9, V ,A), (4) 

where e = l/L and a = tq/t^ are the ratios of the characteristic lengths and times of the activator and the inhibitor, 
respectively. The boundary conditions for Eqs. (Q) and (Q) may be neutral or periodic. 

Kerner and Osipov developed the qualitative theory of large-amplitude patterns in reaction-diffusion systems |3£_ 
(for a comprehensive review on the subject see Refs. [f^]])- They showed that the overall type of patterns is determined 
by the values of e and a, and the form of the nullclinc of Eq. (||), that is the dependence r]{9) implicitly determined 
by equation q{9, rj^A) — for a fixed value of A. For many of systems where patterns may form this nullclinc is N- 
or inverted N-like (Fig. |]). 

According to the general qualitative theory fHH, when e <C 1 and a> 1, only static patterns may form in the 
system; when a< 1 and e > 1 only traveling patterns may form; and when both e <C 1 and a <C 1, all types of 
patterns — static, traveling, and pulsating — may form. 

From the mathematical point of view the fact that 9 is activator means that in some range of the system's parameters 
q' g < 0. In N-systems this condition is satisfied for 9q < 9 < 9' (see Fig. |l|). The fact that r\ is inhibitor means that 
the following conditions hold |||| 

Q; > 0, q' n Q' e < 0, (5) 

and in the whole range of the system's parameters the derivatives Q' g , Q 1 ^, and q'^ do not change signs. 

The systems we are considering have a unique homogeneous state 9 = 9h and 77 = 77^ , where 9h and rjh satisfy 

q(6 h ,ri h ,A)=0, Q(9h,rj h ,A) = 0. (6) 

The homogeneous state is stable for A < Ao, where Aq is the point where 9k = 9q (see Fig. [j]) 

As follows from the qualitative theory [HJg], condition e <C 1 is necessary for the existence of AS and other large- 
amplitude patterns in the considered reaction-diffusion systems. This fact allows to use e as a natural small parameter 
and construct asymptotic solutions by means of the singular perturbation theory ||||[$^,f|l| . Kerner and Osipov have 
shown that as e — > 0, a pattern looks like a collection of "hot" (high values of the activator) and "cold" (low values of 
the activator) regions, separated by the walls whose width is of the order of e |5|J^,^7| -fiof . Thus, in the limit e — * any 
pattern in a (i-dimensional N-system can be described as a (d— l)-dimensional manifold S, corresponding to the walls 
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of the pattern, which separates hot and cold regions fl+ and respectively (Fig. |J). In general, S is a collection 
of an (infinite) number of disconnected orientable submanifolds Si . 

Let us introduce the orthogonal curvilinear coordinates around each submanifolds Si. For a point x let pi be the 
distance from x to Si, and the d — 1-dimensional coordinate the projection to the submanifold pi — const. The 
value of pi is assumed to be positive, if x S f2+ and negative otherwise. In the region of size ~ e around each Si the 
variation of the inhibitor in the direction perpendicular to Si and the variation of the activator along Si are negligible 
compared to the variation of the activator in the direction perpendicular to Si P,|6|]47j| . Therefore, in the vicinity of 
Si stationary Eq. (^) can be approximately written as 

e 2 ^ + e 2 K i (p i ,^ = q(Oi,vlA) ) (7) 
dp; dp t 

where iQ(/?i,£j) is the curvature of Si at a point with the curvilinear coordinates pi and The boundary conditions 
for 9i in Eq. (£]) are 

^(-00) = ^, 0i(+oo) = #3, 0((O) = 0j 2) (8) 
where 6** fc satisfy q(9 l sk , rf a ,A) = 0, with 0^ < < 9\ 3 , and ry* satisfies the consistency condition 

/•+<» / d9\ 2 f e = 3 

e 2 J KiM)[j£j d Pi = q(9,r)l,A)de, (9) 

which follows from Eq. (0), if we multiply it by d9i/dpi and integrate over p.;. Of course, the infinities in Eq. (||) 
actually mean that the boundary conditions for Eqs. (^) should be satisfied sufficiently far from S, namely for \p\ ^> e. 
Note that in the equations above appears only as a parameter. 

When K i — ► the solution of Eqs. ([?]) - (Q) naturally transforms to the one-dimensional sharp distribution (inner 
solution) §,|||]|8|. 

e 2 ^=q(d sh (p i ),T 1a ,A), (10) 
uPi 

where t] s is a constant determined by equation 

" 3 q(8, Vs ,A)d6 = 0, (11) 



and 6 s k are constants satisfying q{9 s k,rj Sl A) — 0. The boundary conditions for Eq. (|10|) are given by Eq. (|8j) with 
?7* = n s , 9\ k = 9 sk . Thus, Eqs. (Q) - (||) describe the sharp distributions of the activator around Si and take into 
account the curvature of Si. 

Far from S the characteristic length of the activator variation is of order one, so for \pi\ 3> e the solution of Eqs. (|^) 
and (^) is approximately give n by the smooth distributions (outer solutions) 9f m (x) and rjf m {x) defined for x S f2±, 



respectively, which satisfy [pU6||3 

= Q(9f m ,vf m ,A), q(9f m , V t m ,A) = (12) 
for x € fi±, respectively, with the boundary conditions 

^ m (^) = ^, Mri = ^, ^0 = ^3, ^ m (^) = ^i (13) 
dpi dpi 



for any Xi € Si. Note that the shape of S itself is determined self-consistently via Eqs. (12) and (|13|). 

According to the singular perturbations theory 0,0], for a; S f2± the asymptotic solution of Eqs. (0) and (ji|) is 
given by the following composition of the sharp and smooth distributions 

9 ± (x) = 9f m (x) + Y^Mx) - 9i h3 ), ^(x) = nf m (x) (14) 

i 

where the sign "plus" goes with #* 3 and the sign "minus" goes with 9 l sl . 
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III. GENERAL METHOD FOR CALCULATING INSTABILITIES 



Let us consider the problem linearized about the static solutions of Eqs. (^) and (|4|) with respect to the fluctuations 

50(3!, t) = S9(x)e iut , 5r](x, t) = 5r)(x)e iut . (1 5) 

The equations describing the fluctuations with the frequency w will become 

iau59 = e 2 A59 - q' e (9(x), r](x))59 - £^(0(x), r]{x))5ri, (16) 

iivSr) = ASrj - Q' e (6(x), r){x))59 - Q' v (9(x),r](x))Sr]. (17) 

As was shown in the previous Section, the solutions 9{x) and r)(x) in the form of a static pattern, around which Eqs. 
(||) and (Q) are linearized, are approximately given by Eqs. ( |l4| ) for sufficiently small e. 

According to the general qualitative theory , the stabilization of a pattern occurs due to the damping effect of 
the inhibitor on the fluctuations of the activator. It was shown that only those fluctuations of the activator that are 
localized in the walls of the pattern and lead to their small displacements are dangerous for the pattern's stability. 
To incorporate this fact into our analysis, let us first consider the fluctuations in the vicinity of <Sj with the fixed 
distribution of the inhibitor. Putting Srj = in Eq. ([lq), we may write 



iauSO = -(H$ + Hi + e 2 S l )S9, (18) 

where 

^ = -e 2 ^ + q'e(e sh (pi),V S ), (19) 



H l e = -e 2 K + q' g (0(x),v(x)) ~ q'e(0s h ( Pl ),Vs), (20) 

and the operator Si is the part of the Laplacian acting on the d — 1-dimensional coordinates , evaluated at Si and 
taken with the minus sign; K is the rest of the Laplacian associated with the curvature of iSj. 

The lowest bound eigenstate of the operator Hg is 59q = d9 s h/dpi, which corresponds to the eigenvalue A = 
|p|.^p7| [I(if] . Indeed, if we differentiate Eq. ( fiol ) with respect to pi, we will get 

^dOsh = a (21) 
dpi 

The function S9q = d9 s \ l jdpi has no nodes; therefore, it corresponds to the lowest eigenstate. The eigenvalues 
corresponding to the excited states are all of order one. This means that the excited states, whose eigenvalues are 
all positive and of order one, are highly damped and arc therefore not important for our analysis of the instabilities 

gifji. 

The operators Hi and e 2 Si in Eq. ( jig ) can be treated as perturbations to the operator Hg. Since the unperturbed 
operator Hg has no dependence on £j, we need to introduce the orthogonal basis of the states corresponding to the 
surface modes on Si. As such a basis, we may choose the eigenfunctions of the operator Si. Then the dangerous 
fluctuations 59 s h are linear combinations of the functions S9u, where 

SOu = Mti)^, Sifa = v u <t>u (22) 
dpi 

and 4>i satisfy 

f ti(ti)MZi)<t 1 ~ 1 Zi = 8u'. (23) 

Of course, for each Si there is its own set of (f>i. We will frequently omit the indices such as i and w, wherever it does 
not lead to ambiguities, in order to simplify the notation. 

Up to now we ignored the reaction of the inhibitor on the fluctuations of the activator. According to the general 
qualitative theory, this reaction can be included into our analysis by means of the perturbation theory MM . The main 
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problem here is to correctly find the response of the inhibitor on the dangerous fluctuations of the activator. The 
formal solution of the problem is of no practical use since one has to diagonalize a complicated operator, nor is the 
expansion in the eigenfunctions since one has to consider the whole spectrum of the problem. A way to do this is to 
use the idea of singular perturbation theory and separate 59 into the localized (sharp) and the delocalized (smooth) 
parts (for more rigorous derivation see Appendix |A|) : 

58 = 58 sh + 89 sm . (24) 

Far from S Eq. ([16]) will become 

iau>56 sm = -{q' g )sm,S9 sm - (cQamSf), (25) 
where the subscript "sm" means that the derivatives are evaluated at the smooth distributions, that is 

(q'e)sm = q' g (0 S m(x),r] sm (x)), (q' v ) sm = q' n (O sm (x), rj sm (x)). (26) 



As follows from the general qualitative theory |5|j6J , for all types of instabilities the condition auj <C 1 is satisfied. For 
this reason we may neglect the left-hand side in Eq. ( p5[ ) , and obtain that 

S6 sm = Sr, (27) 

Let us substitute Eq. (24) into Eqs. ( |l6| ) and (|l7l). Using Eqs. (|2l| ) and (p7|), we can rewrite Eq. ( |l6| ) around S, 
and Eq. (17) in the whole space as 

[iau + e 2 S i +H l )59 sh = -q' ri 5i 1 , (28) 



ilu - A + (Q^) s 



s 



We)* 



5r) = -Q'gSOsh- 



(29) 



Note that we neglected the term Hg59 srn in the left-hand side of Eq. (^8|) since it does not contribute in the first 
order of the perturbation theory. Also note that in writing Eq. (|29] ) we replaced the true distributions 9(x) by the 
smooth distributions 9 sm (x). It is easy to see that this replacement gives negligible difference in 5r]. 
Let us solve Eq. (E^) for Si] by means of the Green's function 



where G(x, x') satisfies 



St 1 (x) = - / Q' 9 (x')G(x,x')69 ah (x')dx', 



[iu - A + C + V(x)} G(x, x') = S(x - x'), 



(30) 



(31) 



where 



and 



C = Q'J9 h ,r lh ) 



q' v (0h,Vh)Qg{9h,Vh) 

q'e( 9 h,Vh) 



V(x) = 



(Q' v )s 



(q' n )sm(Q' e ) 

{Qq ) sm 



c. 



(32) 



(33) 



In view of Eq. (M), the right-hand side of Eq. (]2q) can be regarded as an operator R acting on 59 s h 



R[59 sh ] = q'Jx) / Q' g (x')G(x,x')S9 sh (x')d d x'. 



(34) 



Since the sharp fluctuation 59 s h is the linear combination of the functions 59u defined in Eq. (03) , the integral in Eq. 
(EmI) can be easily calculated. Taking into account that d9 s h/dpi are close to delta- functions PH, we may write 
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R[50 a ] = (Q(O s3 , Vs ) - Q{6 sll i ls )) q 'Jx) / G{x,ii)mi)d d - 1 ^ (35) 



where £j denotes the point on Si and the integration is over <!v The matrix elements of R can be calculated analogously. 
The result of the calculation is 

(i'l'\R\il) = —B f [ G(&,Z i )cf>t,(Z i ,)cf> l (Z i )d d -%d d - 1 Zi ) (36) 

J Si, J Si 

where 

1-6,3 

B = -(Q(e s3 ,v s )-Q(O s i,Vs)) q' v (d,vs)de. (37) 

is a constant depending on A only. Note that in accordance with Eq. (||) the value of B is positive. 

When the distance between the different Si is much greater than e, the overlap between the different 50u is negligible, 
so the operator Hg is diagonal in the z-indices. Then, in the first order of the perturbation theory Eqs. (|2q), (|34|), 
and ( |3li| ) reduce to 

(taLu + e 2 v H )6 u ,S w = eZ-'iii'l'lRlil) - S^iil'lH^tl)], (38) 

where 

Note that the value of Z is of order one since the characteristic length of the activator variation is e. 

Eq. ( p8| ) is the principal equation which determines the stability of an arbitrary d-dimensional pattern in N-systems. 
This equation was derived with an accuracy to e <C 1 and eK max <C 1, where K max is the maximum curvature of a 
given pattern. 

If a pattern possesses certain symmetries, the operators R and Hg are diagonal in the Z-indices. In this case the 
operators in Eq. (|3^ ) can be easily diagonalized (see the following sections). The main problem is to find the Green's 
function G(x,x'). Once this is done, we can obtain the "dispersion relation", which relates uj with the values of A, e 
and a for different types of fluctuations. 



IV. STATIC ONE-DIMENSIONAL AUTOSOLITON IN TWO AND THREE DIMENSIONS 

Let us apply the procedure developed in the previous section to the simplest pattern — static one-dimensional 
autosoliton in two or three dimensions (Fig. ||). Since the AS is localized, the distributions of the activator and the 
inhibitor on its periphery go to the stable homogeneous state — Oh and rj = r]h, where Oh and rjh are determined by 
Eq. (||). In this case, according to Eqs. (||) and (|32|), the value of C > since for Oh < Oq the value of q' g {0h,r]h) > 

For a one-dimensional AS the manifold S consists of two parallel planes where the AS walls are localized. We can 
choose the coordinate directions in such a way that these planes are perpendicular to the z-axis. Then the solution 
for the AS will depend only on z. 

Since the considered static 1-d AS is symmetric with respect to its center Jq,p|, we can assume that the positions 
of its left and its right walls are z\ = —£ s /2 and zi = C s /2, respectively, where C s is the distance between the 
walls. Note that the value of C s itself can be used as a bifurcation parameter instead of A, since there is a one-to-one 
correspondence between them in the whole region of AS existence ^Hj. For this reason, here and further we will use 
the AS width C s as the bifurcation parameter instead of A. 

Since the 6>i,2 are flat, the curvilinear coordinates p\ t 2 coincide with z — z\ and — [z — Z2), respectively (the signs are 
consistent with our definition of pi given in Sec. II), and the coordinates £j coincide with the rest of the coordinates 
of space. Of course, K = 0. 

According to our procedure, let us first look at the operator Si. Here Si x 2 — 2 is the Laplacian, acting on the 

local coordinates on Si^- The eigenfunctions of this operator are just plane waves with the wavevector k along S\ } 2, 
whose eigenvalues are v% — k 2 . 

If we substitute the eigenfunctions of Si t 2 into Eq. ([56|) and use the fact that the system has translational invariance 
in the ^-directions, by integrating over £j-s we will get 
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(l,k\R\l,k') 



BGh 



T' ~2 



S(k-k'), 



r\..l.-\R\2.l/) = -BG, ( "Y> Y ' '* (A '~ //K 



(40) 
(41) 



where Gk(z, z') is the Fourier-transform of a;') in £. Because the AS is symmetric with respect to its center, the 
values of the matrix elements satisfy (1, k\R\l, k) = (2, fc|i?|2, fc) and (1, k\R\2, k) = (2,k\R\l,k). Thus, the matrix 
elements of the operator R in this case can be expressed in terms of the values of the Fourier-transformed Green's 
function at the particular points of the z-axis. 

Let us now turn to the operator Hg. As was shown in the general qualitative theory |^,^|, the functions d9 s h/dpi t 2 
decay exponentially at distances much larger than e. Because of this, for C s ^> e log(e _1 ) the overlap of these functions 
can be neglected. Following the notation of Refs. H§|, we write the diagonal elements of H@ as 

(l,k\H%\l,k') = (2,k\B^\2,k')^e~ 1 ZX S(k-k'), (42) 

where Ao is of order e. The operator in the right-hand side of Eq. (p8| ) can be trivially diagonalized. Since the AS 
possesses central symmetry, the eigenstates of this operator are the symmetric and the antisymmetric combinations 
of 56ik and 502k- These states correspond to the symmetric and the antisymmetric deformations of the AS walls. As 
a result, introducing the functions 



from Eq. 



R 1 (k 1 Lo) = G k 
we obtain the dispersion relation 

„•„., . , ,2 7,2 



C s 


C s ^ 




C s 


2 ' 


2 j 


' 2 


C s 


C s ^ 






2 ' 


2 j 




h A 




eBZ ^Rq i 





(43) 
(44) 

(45) 



Here the subscript "0" corresponds to the symmetric, and the subscript "1" corresponds to the antisymmetric modes 
of fluctuations. The value of Ao can be calculated indirectly. Since the AS possesses translational invariance in the 
z-direction, Eq. (45) should be identically satisfied for k = and to — for the antisymmetric fluctuation. This 
immediately means that 



An 



-eBZ' 1 ^(0,0). 



(46) 



Eqs. ( p3[ ) and (fli| ) define the functions Ra(k, to) and R\{k, u), which describe the inhibitor reaction on the dangerous 
fluctuations of the activator. In general the potential V(z) in Eq. (|3l| ) which determines the Green's function is some 
non-trivial function of z, so Ro,i(k,u)) are some complicated functions of C s . In the case of piecewise-linear model 
one can calculate the values of Ro(k,u>) and R\{k,u) explicitly (see Appendix [b]). One can see that the dispersion 
relations for the fluctuations obtained in this case are identical to those found earlier by Ohta, Mimura, and Kobayashi 
in Ref. |^| who used a different approach. 

To calculate the critical values of A and the parameters of the critical fluctuations in general, we need to know the 
detailed form of the Green's function. According to Eq. (31), the Fourier transformed Green's function Gk{z, z') is 
governed by 



dz- 



+ iu + k 2 + C + V(z) 



G k (z,z')=6(z-z'). 



(47) 



As was said earlier, the potential in the operator in the left-hand side of Eq. ( |47| ) is some complicated function of 
z. However, the problem is greatly simplified for finding instabilities, since, as was shown in the general qualitative 
theory , most of the instabilities occur when C s <C 1 (see also the results below) . This allows to use the value of 
C s as a small parameter and expand the functions Ra t i(k, uS) in terms of it. 

Regarding all this, we are now able to construct the perturbation expansion for the Green's function, considering 
the potential V(z) in Eq. (E?l) as a perturbation. The unperturbed Green's function satisfies 



dz 2 



+ k 2 + iui + C 



G k 0) (z,z')=5(z-z'), 



(48) 



7 



The solution of Eq. (1481) is well known 



(o) expWC+£T£^i]}_ (4g) 

k K ' 2VC + k' 2 + iu V 1 



To calculate the corrections to the Green's function, we use a formula 

p+oo 



/-t-oo 
V{z")Gf{z,z")G^- l \z",z')dz", (50) 
-oo 

where Gj?\z, z') is the ri-th correction. 

Let us denote the contributions from G^™' to Roi(k, u>) as R^i(k,tj), respectively. Since for small values of C s the 
coefficients B, C, and Z only weakly depend on A, we may replace them by their values at A — Ab, where the AS 
size becomes formally zero when e -> §§]. Then the functions R ( °\k,uj) and R ( " ] (k 7 uj) can be written as 



li°\k,u) = — = 1 . {l + cxp(-/Wg + fc 2 + H}, (51) 
2VC- + Ar + «w L J 



^(*,a,) = ^-^±^ + 0(r»). (52) 

Note that R^ — 0(£"), and because of the central symmetry of the potential R± = 0(C 2n+1 ). 
Substituting these values of Ra t i(k, oj) into Eq. (|46|), we will obtain that the leading term of Ao is 

A o = g ■ ^ 53 ) 

Having calculated the value of Ao, we may write the dispersion relation for the symmetric fluctuations 



2 fc 2 = eBZ- 1 - - 1 . {1 + exp(-£ sV / cTF+^)} - R^\k, «)} + 0(/*). (54) 



iaio + e 

As was shown in the general qualitative theory [||,^ 39 4^] , when a is big enough, for C s > C c \ the one-dimensional 
AS becomes unstable with respect to the fluctuation with Re lu — and k = k c 3> 1 , corresponding to the corrugation 
of the AS walls. Note that, according to Eq. (|50|), for k c ~^> 1 the value of R^p (k, u>) which is of order C s , contains a 
small factor oc fc~ 2 , so it can be neglected. Let us calculate the values of L c \ and k c . Putting ui — and neglecting C 
in comparison with k 2 in Eq. (|54|), we will get a transcendent equation, which can be solved for small values of C s . 
The result is 

eZ\- 1/3 . „. /e^ 1/3 



* c = 0.71 1-1 , C c , = 2.64 1-1 . (55) 

Note that the dependences of £ c i and k c on e coincide with those obtained in the qualitative theory for C s <C 1 . 
In the case a <C e there is an instability for £ s > C u with respect to the fluctuation describing the pulsations of the 

AS with the frequency u> = w c >■ 1, and fc = ||| 41 S|. As before, the term R^\k,u) contains the small factor 



and can be neglected. Solving the transcendent equation obtained in this case, we have 

rr\ -2/3 / ry\ 1/3 

aZ \ „ / aZ N 



w c = 0.73l— j , = 0.96 I— . (56) 

Again, the calculated dependences of C u and oj c on the ratio a/e coincide with those obtained in the qualitative 
theory ]^||, and in the analytic studies of the piecewise- linear model |p5| . Also, comparing C c \ and C w , one can see 
that the pulsating instability occurs before the corrugation instability if a < 21e 2 . 

When the size of the AS is greater than the critical size determined by Eq. (pq), the increment of the growing 
fluctuations may be very small. Indeed, according to Eq. (p4[), for C s C c \ we obtain that the increment of the 
most dangerous fluctuations is 7 ~ eB£ s /(2aZ). 
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Now let us turn to the antisymmetric fluctuations. According to Eq. (J52|), the first term in R\{k,ui) which depends 
on lo and k is of order C 2 . As was mentioned earlier, the first correction R\ {k,Lo) is of the order C 3 , so it can be 
neglected. Then Eq. ( (45|) for the antisymmetric fluctuations can be written as 



iaoj + e k 



2;„2 



eBZ^C 2 



\yfC + k 2 + iuu - Vc\ + 0{Cl). (57) 



According to the general qualitative theory |^||, there are two types of antisymmetric instabilities: wriggling of 
the AS walls, and formation of a traveling AS. The first instability is realized when Re lo — 0, k — > and C s > C C 2- 
The second is realized when a <C 1, k = 0, and C s > Ct- We may expand the the right-hand side of Eq. ( p5^ ) in the 
powers of k and to and obtain the following expression 

iaia + e 2 k 2 = ebC 2 s (k 2 +iu>), (58) 

where we retained only the first nonvanishing terms. The constant b is given by 

b = ^m- ^ 

One can see from Eq. (|58|) that the instability Im lo < occurs at k = for C s > Ct, where 
or for i?e lo = and £ s > £ C 2, where 



1/2 



£ c2 = y . (6i) 

Comparing these two formulas, one can see that traveling AS forms before the walls of the static AS become unstable 
with respect to wriggling, when a < e 2 . This fact is also an obvious consequence of an additional symmetry present 
in Eqs. (|) and (|) at a = e 2 . 

Although for C s > C C 2 the AS is unstable, the increment of the most dangerous fluctuations may be extremely 
small. Indeed, according to Eq. (|57|), for £ c2 <C Ct and C C 2 <C C s C c \ we have 

( C 2 B\ 2 C 2 B 
Imax = a' 1 ( j < 1 for k = k max = > 1. (62) 



Comparison of the expressions in Eqs. (55) and ( pl|) for the critical values of C s for symmetric and antisymmetric 
fluctuations shows that for £< 1 the wriggling instability always emerges before the corrugation instability. Likewise, 
according to Eqs. (|56| ) and (|o|), for sufficiently small ratios a/e traveling AS forms before the pulsating one. Notice 
that these general conclusions are in agreement with the analytic investigation of AS in the piecewise-linear model 

EH- 

When the size of the AS becomes comparable with e, the overlap between the eigenfunctions of the operator Hg 
becomes significant, what leads to the instability of the AS. As was shown in the general qualitative theory plpl, the 
instability of an AS with C s ~ e occurs with respect to symmetric fluctuations. Since the asymptotic behavior of 
the sharp solutions is exponential, an extra piece in Eq. ( J54| ) from the operator Hg will have the form aexp(— C s /l), 
where I = eq'g(6 S 3, Vs)' 1 ^ 2 , a is some constant of order 1 ]5||i|.p7| [Io| . Then we obtain that for a ^ e the instability 
occurs with respect to the fluctuations with Re lo — and 

k = k c = 2~ 1 / 3 ( e —\ , C s <£ cb = -loge- 1 , (63) 
whereas for a e the instability is realized with respect to the fluctuations with k — and 

^gj , C s <C bu] = --logae 2 . (64) 

Eqs. d63|) and ( |64] ) were calculated with the logarithmic accuracy and coincide with those obtained in the qualitative 
theory^B. 
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In one-dimensional systems with a ^> e the value of Cb at which the AS disappears will be slightly different from C c b, 
since the only possible value for k there is zero. Putting k = in the dispersion relation we obtain that Cb — Hoge -1 

Up to now we considered the AS whose width C s -C 1. This is justified in 2- or 3-dimensional systems since, 
according to Eqs. ( |55| ) and ([^), the AS becomes unstable when C s <C 1. In one dimension, however, an AS remains 
stable up to the values of C s ~ 1, if a e. In this situation, according to the general qualitative theory HD, there is 
another effect causing the instability. Namely, when C s is close to some Cd ~ 1 at which the activator in the AS center 
reaches the value of 6' , such that q'g(0' o ,T]' o ) = 0, the solution in the form of an AS disappears. As a result, a local 
breakdown occurs at the AS center, causing the AS to split, so that eventually the whole system becomes filled with 
a complex pattern ||||. This can also be seen from Eq. ([45]), if we take into account that at the point z = where 
= 0' the potential V(z) in Eq. (^) becomes singular. In the case a < e the values of £ w and Ct corresponding 
to the instabilities leading to formation of pulsating and traveling AS, respectively, may be less than Cd, so the AS 
cannot reach the point where the local breakdown occurs. There is also a possibility that the point A = Ad where 
C s = Cd is preceded by the point A = A c where the homogeneous state of the system becomes unstable |||| . In this 
case the periphery of the AS becomes unstable. This can also be seen from Eq. (^) if we take into account that, 
according to Eq. (^), for A > Aq the value of C becomes negative and therefore the tails of the Green's function 
will become oscillatory. When a gets larger, the right-hand side of Eq. ( [45| ) always remains of order one, so at some 
critical values of a > e the instability leading to the formation of traveling and pulsating AS disappears. 



V. HIGHER-DIMENSIONAL RADIALLY-SYMMETRIC AUTOSOLITONS 



Let us now turn to higher-dimensional AS. First we consider spherically symmetric AS in three dimensions. In this 
case there is only one manifold S which is simply a sphere of radius 1Z S . As a set of orthogonal curvilinear coordinates 
we choose the usual spherical coordinates p, d and tp, except p will be measured from the surface of the sphere rather 
than from the origin, in order to be consistent with our initial definition. 

The operator S in the the considered case is 

S = j- 1 ^ sin^A + * |U . (65) 

|_ sin v ov ov sin -d oip z J 

The eigenfunctions of the operator S are just spherical harmonics cj>i m = Y; m (i?, tp)/lZ s with the eigenvalues v\ = 
1(1 + 1)/TZ^. The factor 1/1Z S in the eigenfunction ensures the proper normalization. 

Now let us calculate the matrix elements. First of all, since the system possesses spherical symmetry, the only 
nonvanishing matrix elements of the operator H@ are the diagonal elements, which are all equal to each other. Due 
to the same symmetry, the operator R is also diagonal, and its diagonal elements are independent of m. According 
to Eq. (|36|), the diagonal matrix elements of R are 

(lm\R\lm) = -BRi(lu), (66) 

where 

R 1 (uj)= f [G(n s ,tf,p;n s ,tf\<p')YC m (#,cp)YU#\v')K 2 s do'do, (67) 



do is the element of the solid angle, and the Green's function is written in terms of the spherical coordinates. 
To calculate Ri(lo) let us note that Ri(u>) = 1Z^Gi(1Z s ,1Z s ), where 

Gi(ry)= f I ' G{r,d^-r'J>,p>)YC m {d,<p)Y lm {{)',pi)dodo' (68) 



is the Green's function satisfying 



r 2 



d 2 2d 1(1 + 1) 
dr 2 r dr r 2 



+ C + iu + V(r) 



Gi(r,r') = S(r-r'), (69) 



Eq. ( |69| ) follows from Eq. (|3lj) if we first rewrite it in the spherical coordinates r, ip, multiply both sides by 
<p)Yim('&' i V 3 ') an d then integrate over do and do' , taking into account the hermiticity of the angular part of the 
Laplacian and the orthonormality of the spherical harmonics. 
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Using the above mentioned properties, we can now write Eq. (|38|) for the spherically-symmetric AS of radius 7Z S 
in the form 

iau + 6 ^ ^ + A = -eBZ- l Ri{uj), (70) 

where Ao is the contribution from Hg . Eq. (|7^) is the dispersion relation for a fluctuation characterized by the number 
I. 

It can be shown by direct calculation that to the first order in e/lZ s the value of Ao is zero. This means that in order 
to calculate Ao from its definition one should know the distributions 9{x) and rj{x) for the AS with greater accuracy 
than that of Eq. (fu]). Also, the second order of the perturbation theory has to be taken into account. However, 
these difficulties can be avoided, if we use the fact that the system possesses translational invariance. Indeed, Eq. 
( |67| ) should be identically satisfied for I = 1 and v = 0, so we immediately obtain that 

^o = -^-eBZ- 1 R 1 (0). (71) 

As we will show below, the instability of a spherically symmetric AS occurs at 1Z S <C 1. To find the function Ri(u>) 
let us use the idea of the previous Section and seek for the Green's function of Eq. ( |69[ ) , treating the potential V(r) 
as a perturbation. The zeroth order Green's function is well known 

. Il+l/2(Kr')K l + 1/2 (nr") , ^„ 

G i (r',r")={ i l+1/2 ( K ^)^ l+1/2 ( K r') r , >r „ ( 72 ) 

y/r'r" ' ' 

where n = y/C + icu; Ki+U2(z) an d h+i/2( z ) are the modified Bessel functions. The corrections to G*f\r' , r") are 
then given by 

POO 

G\ n) (r', r") = - / V(r)Gf ] (/, r)G ; ( "" 1) (r, r")r 2 dr. (73) 
Jo 

As a result, from Eq. d^) we obtain that 

Ri°\uj) = 1I s I i+1/2 (k1I s )K i+1/2 (k1I s ), (74) 

and, as can be seen from Eq. ©, R[ n) (cu) = ll 2 s G<l n) (1l s , K s ) = 0{K 2 s n+1 ). 

As before, we will use the values of B, C, and Z evaluated at the point A — where 1Z S — > as e — > 0. Expanding 
the Bessel functions at 7?. s < 1, we may write Eq. ( frc| ) as 

iao; + £ 2 K; 2 (! + 2)(i - 1) = eS^- 1 ^ Q - . (75) 

Eq. ( [75| ) describes the fluctuations in the case |w|7?.^ <C 1. This is satisfied for the thresholds of the aperiodic 
instabilities. Simple calculation shows that an AS becomes unstable with respect to the aperiodic I = mode when 

K s < K b = {^pj 1 . (76) 

For 1Z S > lZt> it becomes unstable with respect to the aperiodic fluctuations with I > 1 at 

/3(Z + 2)(2; + l)eZ\ 1/3 , . 

Us > Ucl = { 2B ) ■ (77) 

One can see from this equation that the first instability point corresponds to I = 2: 

^ c2 = (^-g-J • (78) 

Thus, a spherically-symmetric AS can be stable only if its radius satisfies IZb < Ti-s < ^c2, where IZb and 1Z C 2 are 
given by Eqs. ([76]) and (]7q), respectively. 
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The I = 1 mode corresponds to the translation of the AS as a whole. As in the case of the one-dimensional AS, for 
some value of a <C 1 the spherically-symmetric AS becomes unstable with respect to the fluctuation leading to the 
formation of the traveling AS. The instability with respect to the I = 1 mode occurs for some TZ S > TZt when u> — > 0. 
Expanding the Bessel functions in the right-hand side of Eq. d7G) for small TZ S and u>, we get 



2 TZ 

iauj = eBZ- l n A s { ioj— - (ioj) 2 '-== \ + h. o. terms. (79) 



15 24VC 

According to this equation, the static AS transforms into the traveling one, when TZ S > TZt, where 

/15aZ\ 1/3 

TZt = -5-5- ■ (80) 



V 2eB 

Now let us study the instabilities with Re lo ^ 0. Consider an AS stable when a>l, and whose radius is therefore 
of order e 1 / 3 . As follows from Eq. (|70"|), in order for an instability to occur the frequency Be u> at the threshold of the 
instability should be big enough, so that the argument of the Bessel functions in Eq. (|74|) is of order one. This means 
that lo ~ TZj 2 ~ e~ 2 / 3 and therefore the critical values of a are of order e 2 . Indeed, let us introduce the new variables 

eZ\ 2/3 BTZl 



a = a/e\ Co = u ^— J , P = ~rf ■ (81) 

Substituting these variables into Eqs. (70), ( |7l| ) and (|74|), after some algebra we obtain the following transcendent 
equation which has explicit dependence on a, p, and ui only. 

lauo + (l + 2)(l- 1) P -^ = pV3 Q _ h+1/2 {pi/s^i/^V?} K l+1/2 {pVWVif) . (82) 
We solved Eq. (82) numerically for I = and 2 in the region of p where the AS is stable with respect to the aperiodic 



fluctuations. The resulting stability diagram is shown in Fig. (^). The rescaled critical frequency as a function of 
a for I — is also presented in Fig. (||). From Fig. (Q) it is clear that when a gets smaller, the AS always loses 
stability with respect to the I — pulsations first. The AS is always unstable if a < a c — 4.4e 2 . For a > 6.7e 2 the AS 
destabilizes with respect to the aperiodic / = 2 mode first, if its radius is increased. All other instabilities, including 
the one leading to the formation of a traveling AS, occur at smaller values of a and are, therefore, secondary. 

Let us now consider the case of radially-symmetric AS in two and three dimensions. Because of the close analogy 
with the case of spherically-symmetric AS, we will only outline the derivations, focusing mainly on the obtained 
instability criteria. 

If r, if and z are the cylindrical coordinates, then the coordinate p = — (r — TZ S ), and the operator S is 

1 d 2 d 2 

S = ~n 2 d^~d^- (83) 



The eigenfunctions of this operator are 4>km — (47r 2 7?. s ) 1 / 2 e lkz + lmi P with the eigenvalues Vkm 
The non- vanishing matrix elements of the response operator in this case are 



iz 2 



(mk\R\mk') = -BB m (k, oj)5(k - k'), (84) 
where the function B m (k,uj) can be expressed in terms of the Green's function given by 

d 2 1 d 



—To - - J f- + 1 ^r + C" + /s 2 + w + V(r) 
ar r ar r z 



G km (r,r')=S(r-r'), (85) 



as R m {k,uj) = TZ s G km {TZ s ,TZ s ). 

The dispersion relation for the fluctuations with particular values of k and m, obtained from Eq. (|38|), is 

2 2 

iauj + e 2 k 2 + + A = -eBZ^Rmik, lu). (86) 

Because of the translational invariance, this equation should be satisfied identically for k — 0, lo — and m = 1. This 
gives us 
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e 2 

X = --^-eBZ- 1 R 1 (0,0). (87) 

As before, for small values of 1Z S we will seek for the function R m (k,oj) perturbatively. The zeroth order Green's 
function here is 



(J J>\ - / ImW)K m (nr"), r' < r", 
Gfcml ' ' ~ \ I m W")I m (Kr'), r' > r", 



where k — y/C + k 2 + icu, I m (z) and K m (z) are the modified Bessel functions. The corrections to the Green's function 
are given by 

/>oo 

G { :l(r'y') = - / V^G^r'^Gt^iryydr. (89) 



From this we find that 



R$(k,uj) = K s I m (KK s )K m ( K K s ), (90) 



and that the corrections to R m {k,uj) from are R^\k,uj) = n a G ( ^l(Tl s ,Tl s ) = 0(TZ 2 s n+1 ). 

Let us study the stability of the radially-symmetric AS in two dimensions first. Putting k = in Eq. (|S6| ) with 
R m (k,uj) given by Eq. (pp|), we obtain that for m > 1 and 7?. s <C 1 the aperiodic instability occurs at 1Z S > lZ cm , 
where 

K„ = ( 2em(m B + " z )' A (91) 

According to this equation, when the value of 1Z S is increased the AS becomes unstable with respect to the fluctuation 
with m = 2 when 1Z S > TZ C 2, where 

12eZ^ 1/3 



n c2 = I — j . (92) 

When the value of 1Z S is decreased, at 1Z S < IZj, the radially-symmetric AS in two dimensions becomes unstable 
with respect to the aperiodic fluctuations with m = 0. According to Eq. (p^), for small values of 1Z S we have 

eZ N 1/3 



^=^j - (93) 

where 6i = — log(1.477?. s V / C'). Since 6i only weakly depends on 1Z S , for 1Z S ~ 0.1 we may put &i ~ 2, and for 1Z S ~ 0.01 
we may put &i ~ 4. Thus, a static radially-symmetric AS in two dimensions can be stable only if TZb < TZ S < T^c2 
with TZb and 1Z C 2 given in Eqs. (|9^) and (^2|), respectively. 

When a gets sufficiently small, the radially-symmetric AS in two dimensions typically destabilizes with respect to 
the m = pulsations first. If we introduce the variables of Eq. ( |3l"l ) into Eq. (86), for small 1Z S we will obtain an 



equation similar to Eq. (p2|), which depends on p, a, and ui only. Solving this equation numerically for m = and 
2, (the case m = 1 will be discussed below) we obtain the stability diagram for the radially-symmetric AS in two 
dimensions (Fig. ||). Figure || shows that when a < a c = 6.1e 2 , the AS is unstable for all values of TZ S . When 
a > lie 2 the AS first destabilizes with respect to the aperiodic fluctuations with rn = 2. The dependence of the 
rescaled critical frequency ui on a for the m = pulsations is presented in Fig. ^. All instabilities with m > 2 occur 
at smaller values of a, so they are secondary. 

Now let us turn to the radially-symmetric AS in three dimensions. Solving the dispersion equation for small 1Z S for 
m = 0, we obtain that the AS becomes unstable when 1Z S <lZb, where 

/2.25eZ^ 1/3 
Kb = ( 

or when 1Z S > 1Z C , where 



,7.5eZ\ 1/3 

^c=(-g-J , (95) 
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when 



4.1e^ - 1/3 



k = k c = I — j . (96) 

Comparing Eq. (|9^) with (^2|), one can see that in contrast to the radially-symmetric AS in two dimensions, when 
1Z S is increased, the radially-symmetric AS in three dimensions destabilizes with respect to m — mode first. 

In the case m = 1 an AS destabilizes with respect to the fluctuations with small k at Re uj = 0. Expanding the 
Bessel functions in the right-hand side of Eq. ( J86| ) , we will obtain the following equation 

iauo + e 2 k 2 = eBZ- 1 b 2 'Rl{k 2 +iuj), (97) 

where b 2 — — log(1.147?. s v / C')/4. The value of b 2 weakly depends on TZ S , so for 1Z S ~ 0.1 we may put b 2 — 0.5, and 
for 1Z S ~ 0.01 we may put 62 — 1- It can be seen from Eq. (|97j ) that when a » e 2 an AS becomes unstable when 
7Z S > lZ c i at k — > 0, where 

eZ N 1/3 



**=UfeJ ' (98) 

Comparing Eq. d98| ) with Eq. (63) one can see that for e < 10~ 3 we have 1Z c i < TZb, and, therefore, the cylindrically- 
symmetric AS is always unstable. However, note that the increment of the fluctuations with m = 1 and small k may 
be extremely small. 

Similarly, as follows from Eq. (|97|), when a ~ e 2 the AS becomes unstable for 7£ s > at k = 0, where 



and transforms into the traveling AS. As we already noticed, when TZ S is not very small, the instability leading to the 
formation of a traveling AS occurs at smaller values of a than the instability with respect to the m — pulsations. 
However, when 1Z S < 0.01 the coefficient b 2 in Eq. ( |99"| ) becomes such that the line p(a) in Fig. @ corresponding to 
the instability with m — 1 crosses the curve corresponding to the threshold of the m = pulsations when p < 12. 
In this situation a radially-symmetric AS in two or three dimensions may become unstable and transform into the 
traveling AS before it destabilizes with respect to the m — pulsations. However, this is a rather unrealistic situation 
since in order for an AS to have a radius 1Z S < 0.01 and be unstable with respect to the m = 1 mode, one should 
10~ 6 and a < 10~ 12 . 



VI. CONCLUSION 



Thus, in the present paper we developed an asymptotic theory of the instabilities of arbitrary d-dimensional static 
patterns which may form in a wide class of reaction-diffusion systems of the activator-inhibitor type. This theory is 
based on the natural smallness of the parameter e = l/L, where I and L are the characteristic length scales of the 
activator 9 and the inhibitor 77, respectively. In fact, as was already mentioned, if the length scales I and L, as well 
as the normalization of 9 and rj, are chosen properly, the condition e <C 1 is not only sufficient, but also necessary for 
the existence of AS and other large amplitude patterns |^,^]. 

Within the presented theory we analyzed different types of spontaneous transformations of the simplest static 
patterns in two- and three-dimensional systems into more complex static, pulsating, and traveling patterns. We 
showed that the criteria corresponding to these transformations are universal in the sense that they are practically 
independent of the specific nonlinearities of the system and are determined only by the two parameters: e and a; and 
three numerical constants: B, C, and Z, which have all necessary information about the nonlinearities. If the length 
scales and the normalization of 9 and 77 are chosen properly, the constants B, C, and Z are necessarily of order one, 
and the constant B can in principle be small. 

Let us summarize the results of our analysis. 

According to the formulas of Sec. IV, when a 3> e a one-dimensional AS in two and three dimensions (stripe) is 
less stable than the AS in one dimension. As follows from Eqs. ( p5| ) and (|6l|), for e < 1 we have C c2 < C c \. This 
means that when the control parameter A is increased and the AS widens, the AS always destabilizes with respect 
to the wriggling of its walls first. It is natural to expect that for C c2 < C s < C c \ a one-dimensional AS will deform 
into a twisted band that fills the whole volume of the system. When A is further increased so that C s > £ c i, each 
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wall of the pattern becomes unstable with respect to the fluctuations with the characteristic wave vectors k c ~ e -1 / 3 . 
As a result of the development of this instability, fingers will start to grow from the walls of the pattern. Eventually, 
the volume of the system will become filled with a labyrinthine pattern. The instability will persist until the distance 
between the walls and their curvature radii become of order e 1 / 3 . This phenomenon was observed recently in the 
experiments by Lee and Swinney |nj and in the numerical simulations of a two-dimensional reaction-diffusion system 

When A is decreased and the AS narrows, at C s < C c b it destabilizes with respect to the corrugation of its walls 
with the wavevector k c ~ e^ 1 / 3 . As a result, granulation of the static AS will occur and eventually the resulting 
granules with the radius of order e will disappear. 

According to Eqs. ( |6l| ) and (|60|), for a < e 2 we have Ct < C c %- This means that as C s is increased, a one- 
dimensional AS will always transform to a traveling band first. Note that the condition a < e 2 here is exact and does 
not depend in any way on the nonlinearities of the system. 

As follows from Eqs. ( p6| ) and (p0|), there may be a rather wide range of the parameters a and e for which the 
pulsation instability emerges before the instability leading to the formation of traveling AS as the width of the AS is 
increased or the value of a is decreased. Indeed, the condition C u < Ct is satisfied when 



where 



a > ep u , (100) 



A- = 65oi^' (101) 



Note the huge numerical factor in the denominator of Eq. (101). Because of it the instability with respect to pulsations 



will emerge before the instability to traveling AS in the majority of the real systems. At the same time, as follows 
from Eqs. ( j56| ) and (|6l|), the instability with respect to pulsations is the first, i. e., £ w < C C 2, if 

a < e 5/2 /3j 1/2 . (102) 

These two conditions can be satisfied at the same time only if 

e > [3 U . (103) 

As a result of the instability with respect to pulsations a static AS may transform into a stationary breathing pattern 
(pulsating AS), if the parameters of the system are finely adjusted Js],^). However, as we see from our numerical 
simulations, in most cases the walls of the AS go so far apart that a local breakdown occurs at the AS center, and 
eventually the AS produces two one-dimensional AS traveling in the opposite directions. 

When the width of the AS is decreased, at a < e 2 the AS destabilizes with respect to pulsations and disappears, if 
C s < Cbui, where Cbu> is given in Eq. (M). When the value of a decreases, at a — a c we have Cbuj — Ct- For smaller 



values of a the AS is always unstable. According to Eqs. (60) and (|64|), the value of a c is given by 

a c ~e 3 (loge- 1 ) 2 . (104) 

The transformations and the evolution of the static spherically-symmetric AS in three dimensions and static radially- 
symmetric AS in two dimensions are very similar. These AS are stable only in the relatively narrow range of their 
radii. When a is big enough, the AS is stable if IZb < 1Z S < TZc2, where IZb and 7Z C 2 are of order e 1 / 3 and are given 
by Eqs. ([76]) and (f78|) for the spherically-symmetric, and by Eqs. (|94|) and (|9^) for the radially-symmetric AS in 
two dimensions, respectively. As the control parameter A is decreased, at TZ S < TZb the AS will abruptly disappear. 
As A is increased and the AS radius becomes greater than 7l C 2, the AS loses stability with respect to the radially- 
nonsymmetric fluctuations with I = 2 first. The growth of these fluctuations may lead to the splitting of the AS into 
two, or to the growth of a pattern with the sophisticated geometry. This complex pattern may further get more and 
more complicated as a result of the fingering instability, if the distance between the pattern's walls exceed the value 
of order e 1 / 3 . The process of splitting and complicating will go on until the whole system becomes filled with a very 
sophisticated labyrinthine pattern (Fig. |^) . The pattern should not necessarily be connected because of the possibility 
of splitting. Thus, there is a remarkable phenomenon characteristic to the considered class of nonlinear systems: as a 
result of the instability of the AS localized in a small portion of an extended system the whole system becomes filled with 
a complicated pattern. These conclusions explain the effects of splitting, self-replication, and formation of labyrinthine 
patterns found recently in the experimental and numerical investigations of some two-dimensional reaction-diffusion 
systems ©HH. 
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Static cylindrically-symmetric AS in three dimensions can only be stable if e > 10~ 3 . If this condition is satisfied, 
the cylindrically-symmetric AS will destabilize with respect to wriggling (m = 1 mode) when 1Z S > 1Z C \, where 1Z c i 
is given by Eq. (98) as its radius is increased. If the radius of the AS is decreased, at 1Z S < IZb, where lib is given by 
Eq. (94), the AS will destabilize with respect to the corrugation of its walls (m = and k = k c , where k c is given by 
Eq. (96)). As a result, the AS will granulate and transform into a number of spherically-symmetric AS. 

When the value of a is decreased, a stable radially-symmetric AS in two dimensions and spherically-symmetric 
AS in three dimensions lose their stability with respect to the radially-symmetric fluctuations oscillating with some 
characteristic frequency (see Figs. § - 0). Only the radially-symmetric AS in two dimensions whose radius 7l s < 0.01 
can spontaneously transform to a traveling AS before it destabilizes with respect to the m = pulsations. However, 
as we already mentioned, this situation is possible only for unrealistically small values of a and e, therefore, this 
bifurcation, which was recently discussed in Ref. [^4j is secondary in most of the real reaction-diffusion systems. As a 
result of the instability with respect to the radially-symmetric pulsations the AS may disappear, or, if the parameters 
of the system are finely adjusted, a stationary pulsating radially-symmetric AS may form [^|^]. However, as we see 
from our numerical simulations, in most cases the growth of the amplitude of the AS pulsations leads to the local 
breakdown in the AS center and the formation of a traveling wave in the form of a ring with the radius monotoncly 
growing with time. 

Static radially-symmetric AS of any radius are always unstable if a < 4.4e 2 in three dimensions, or if a < 6.1e 2 in 
two dimensions (in the case of extremely small e and a this value may be greater: see the discussion above). In this 
situation only traveling waves and pulsating patterns will form in the system. Note that these conclusions are totally 
independent of the specific nonlinearities of the system. 

In our analysis we considered only monostable systems. However, the results obtained by us remain true in 
bistable systems as well. In particular, it can be easily seen that a static one-dimensional front connecting two stable 
homogeneous states is always unstable with respect to the fluctuations with the wavevector k ~ e^ 1 / 3 . 

So, we see that our results are universal and applicable to the wide class of physical, chemical, and biological systems 
which can be described by Eqs. ([!]) and (||), if the nullcline of Eq. ([j]) is N- or inverted N-like (Fig. [l]). In such 
N-systems the universality of the obtained results is related to both the form of the nullcline and the smallness of the 
parameter e. At the same time there are systems which are described by Eqs. (^) and (0), in which the nullcline of 
Eq. ([!]) is V- or A- like. For example, the biological morphogenesis system of Gierer and Meinhardt [jllj, the models of 
axiomatic chemical reactions — the Brusselator [Q and the Gray-Scott model |2^] all have V- or A-like nullclines. In 
such V-systems at e <C 1 spike patterns of giant amplitude may form . The properties of such patterns essentially 
differ from those forming in the N-systems we considered here. For this reason it would be incorrect to use the results 
of the present paper to interpret the results of the numerical simulations in V-systems. However, the simulations 
performed by Pearson in the two-dimensional Gray-Scott model showed |2^] that the effects of the granulation of 
static one-dimensional AS, splitting and self-replication leading to the formation of complex patterns which fill the 
whole space are seen in V-systems as well. 



APPENDIX A: THE EIGENVALUE PROBLEM 



Let us consider the exact stability problem for a static pattern. Equations (|16|) and (|17]) can be written in operator 
form as 



where 



and 



H^e = UiSr], 
H 2 Sr] = U 2 S9, 

Hi = iauj — e 2 A + q' g , 
H 2 =iw-A + Q' 



Ul = U 2 = -Q'g. 



Substituting Eq. QA2[) into Eq. flAl| ), we obtain 

(Hi - u 1 H 2 1 u 2 )se = xse, 



(Al) 
(A2) 



(A3) 
(A4) 



(A5) 



(A6) 
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where A should be put to zero. Solving this eigenvalue problem and then requiring that A = 0, we may obtain the 
value of to. In fact, this allows to think of A as an infinitesimally small quantity. 

In the considered problem the operator UiH^V-i may be treated as a perturbation to the operator Hi HI). We 
would like to find the solution of Eq. (AC) corresponding to the lowest eigenvalue. In view of the discussion in 
Sec. Ill, to the zeroth order in e the eigenfunctions of the operator H\ are linear combinations of the functions 

(0) 



59 



/e/Z 59a where 59u are defined in Eq. ( |22| ) and Z is defined in Eq. ( p9[ ) (the coefficient in front of 
ensures the proper normalization), and their corresponding eigenvalues are of order e. It can be easily seen that 



{i'l'\U 1 H^ 1 U 2 \il) ~ e, 



(A7) 



where the matrix element is calculated with the functions 59^\ However, one should be careful in calculating A since 

the matrix element from the bound state 59^ to the state of the continuous spectrum of the operator Hi with the 
wavevector k ~ 1 has the following estimate 



{k\UiH^U 2 \il) 



(A8) 



So, the second and higher order corrections of the perturbation theory given by the transitions from the bound states 
to the long-wave continuous spectrum will be of the same order as the first-order contribution from the diagonal 
element. 



According to Eq. ( A6) and the fact that the unperturbed eigenvalues of the operator Hi are of order e, the improved 
function <50-, which contains the corrections of order ^fl can be written as 



{1 + H^U X H^U 2 + (H^UiH- 1 ^) 2 + . . .}59f t 



(0) 



(A9) 



Of course, as should be in the perturbation theory, the operator H 1 1 actually projects out the 59 a components. If 



we now substitute this function for 59 into Eq. (A6), multiply it by 50^? , and then integrate over the volume of the 
system, to the first order in e we will arrive at the following equation 



(i'l'\Hi\il) - (i'l'\R\il) =XS W 6 W , 



(A10) 



where 



R = UiH^U 2 {l 



H7 l UiH~ 1 U 2 



{H^UiH^U 2 ) 2 + ...}. 



(All) 



One should not confuse the matrix elements in Eq. ( AlOj) with those of Eq. ( |36| ) and ([58]), since they are calculated 
with the eigenfunctions which have a different normalization. 

To the first order in e, we may replace the true distributions of the activator and the inhibitor in the operator R by 

59 u has the characteristic length scale 



59 



the smooth distributions. According to Eq. ( A9), the function 59 s 
1. If we neglect the term au; in Eq. (A3), the operator H^ 1 reduces to ^(^sm^ilsm^)) 



(i) 



Then the definition 



of the operator R in Eq. (All), together with the above mentioned property of the operator Hi, means that in the 
calculation of the inhibitor resp onse one should consider the fluctuations 5rj and 59 sm to be related by Eq. (p5|). With 
all these approximations, Eq. ( |A10| ) is equivalent to Eq. (|38|). 



APPENDIX B: PIECEWISE-LINEAR MODEL 



It seems that the only model in which it is possible to find the exact Green's function of Eq. (|3l]) is the well-known 
piecewise-linear model of a reaction-diffusions system, which is described by the following equations |l5| ] 

a— = e 2 A6> - 6 — rj + H(9 - A), (Bl) 



^L=Ar 1 + 9- 1 r ] , (B2) 



where H(x) is the Heaviside function. 
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The homogeneous state of this system is 9h = 0, r\h = 0. It can be easily verified that the values of the parameters 
describing the sharp distribution of the activator are 8 s i = A — |, 9 s2 = A, 6^3 = A + i, and r/ s = | — A [|l3| , p~6| ~p~8|| . 
In view of Eqs. ([Bl]), (|B2|), @, and @, we obtain that in this model 5 = 1 and C = 1 + 7. 

In order to find Z we need to know the sharp distribution. According to Eq. ([Io|), for this model 

^"U+i-W-p/e), p>0. (B3) 



From Eqs. (B3) and ( p9[ ) we obtain that Z = 1/4. Note that B, C, and Z are just constants independent of A. 

Having calculated the constants B, C, and Z, we can easily find all instability points. Let us consider the one- 
dimensional AS, for example. According to Eq. (55), for a S> e the AS becomes unstable with respect to the 
symmetric fluctuations with 

k c = 1.12e -1/a , at C c = 1.66e 1/3 , (B4) 
whereas, according to Eq. (|56|), in the case a <C e it destabilizes with respect to the fluctuations with 

'2/3 . /Q\l/3 



(jj , at £ w = 0.60 (^J . (B5) 



Eq. ( p35| ) improves the accuracy of the results obtained in Ref. |25|] . 

Similarly, from Eq. ( |59|) one concludes that in this system b = l/(2i/l + 7). According to Eq. (|6l|), in the case 
a>fa one-dimensional AS becomes unstable with respect to the antisymmetric fluctuations with k — > at 

£ s >£ c2 = (l+ 7 ) 1/4 (2e) 1/2 , (B6) 



what agrees with the result of Ref. J47j. In the case a<e, according to Eq. (|6C|), a static AS spontaneously transforms 
into traveling when 

£ S >£ T = (1 + 7) 1/4 (^) 1/2 . (B7) 

This result coincides with the one obtained in Ref. |^5j . 

When C s becomes comparable with e, an AS becomes unstable with respect to the symmetric fluctuations. Ac- 
cording to Eq. (|6^), for a 3> e the instability occurs at C s < Cb with respect to the fluctuations with k = k c , 
where 

k c = 1.266- 1 / 3 , £ b = jloge-\ (B8) 



and, for a <C e, according to Eq. (|64|), at C s < Cbu with respect to the fluctuations with uj — lu c , where 

w c = 2(~J , C bu) = --logae 2 , (B9) 



since in this case I — e. Similarly, in one dimension £5 = eloge -1 . 

As we noted in Sec. IV, there is a one-to-one correspondence between the AS width C s and the control parameter 
A . In this system this correspondence is given by 

A = -LLI — (BIO) 

2(1 + 7) 1 ' 

where 2(1+7) ^ ^ ^ 5' Thus, knowing the critical values of C s , we can easily calculate the values of A at which the 
instabilities occur. 

Finally, the value of /3 U , which appears in Eq. (101) and determines the region where the AS becomes unstable 
with respect to pulsations before it transforms into a traveling AS for this model is 

U =6.2- 10- 3 (1 + 7)" 3/2 - (Bll) 

This formula also improves the accuracy of the results obtained earlier in Ref. |25) . 
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Because of the singular character of the nonlinearity, the potential defined in Eq. ( |33| ) is identically zero, so the 
dispersion relation for the one-dimensional AS of an arbitrary size in this model becomes 

iau; + e 2 k 2 + \o = 26 (l ±e ~Wi+7+fc 2 +^\ ( B i2) 

y/1 + 7 + fc 2 + iuJ I > 

where the sign "+" goes with symmetric fluctuations, whereas the sign "-" goes with the antisymmetric ones. Ac- 
cording to Eq. ([i"6|), 

One can see that Eq. ( |B12| ) coincides with Eq. (5.4b) of Ref. jl7| by a different method. 



As in the case of one-dimensional AS, the unperturbed Green's functions from Eqs. (72) and (|8§|) for the spherically- 
and cylindrically-symmctric AS, respectively, are the exact Green's functions. For this reason the exact dispersion 
relations for the radially-symmetric AS in this model are given by Eqs. ([70]) and (|6|) with the functions Ri(uj) and 
R m (k,oj) given by Eqs. ( pT4| ) and (|90|), respectively. It is easy to see that the dispersion relations obtained in this 
fashion coincide with those obtained in Ref. as well. 
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FIG. 1. 

Qualitative form of the nullclines of Eqs. (0) and 



FIG. 2. 

"Hot" and "cold" regions forming a pattern. The walls of the pattern are localized in the region of order e around 



FIG. 3. 

Distributions of 9 and 77 in the form of a one-dimensional AS. On the AS periphery 9 goes to Oh and 77 goes to r/h- 
Dashed lines indicate the regions where the AS walls are localized. 



FIG. 4. 

Stability diagram for the spherically-symmetric AS in the rescaled variables a and p. The solid curves correspond to 
the thresholds of the instabilities for I = 0, I = 1, and 1 = 2. The bottom and top horizontal lines correspond to the 
thresholds of the aperiodic instabilities for I = and 1 = 2, respectively. 



FIG. 5. 

The rescaled frequency ui vs. a at the threshold of the / = instability of the spherically-symmetric AS. 



FIG. 6. 

Stability diagram for the radially-symmetric AS in two dimensions in the rescaled variables a and p. The solid 
curves correspond to the thresholds of the instabilities for m = 0, m = 1, and m = 2. The top horizontal line 
corresponds to the thresholds of the aperiodic instability for 1 = 2. The bottom horizontal line shows schematically 

the threshold of the aperiodic instability for m = 0. 



FIG. 7. 

The rescaled frequency ui vs. a at the threshold of the m = instability of the radially-symmetric AS in two 

dimensions. 
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FIG. 8. 

Formation of a labyrinthine pattern. Numerical solution of Eqs. (S) and ([IJ) with q — 8 3 — 9 — 77 and Q = 6> + ?/ — /I 
with e = 0.05, a = 0.1, A — —0.3 at different times. The system size is 20 x 20. At t — the homogeneous state of 
the system is excited in the region 0.5 x 0.7. Distributions of the activator at times t — 0, 16.5, 30, 65, 100, 300. 
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